MOLECULAR ECOLOGY 



Molecular Ecology (2013) 22, 3476-3494 



doi: 10.1111/mec.l2341 



Combined analyses of kinship and F ST suggest potential 
drivers of chaotic genetic patchiness in high gene-flow 
populations 



MATTHEW IACCHEI/t TAL BEN-HORIN, J KIMBERLY A. SELKOE,*§ CHRISTOPHER E. BIRD/K 
FRANCISCO J. GARCIA- RODRIGUEZ** and ROBERT J. TOONEN* 

*Hauoai'i Institute of Marine Biology, School of Ocean and Earth Science and Technology , University of Hawai'i at Manoa, 
Kdne'ohe, HI 96744, USA, ^Department of Biology, University of Hawai'i at Manoa, Honolulu, HI 96822, USA, %Bren School 
of Environmental Science and Management, University of California, Santa Barbara, Santa Barbara, CA 93106, USA, ^National 
Center for Ecological Analysis and Synthesis, Santa Barbara, CA 93101, USA, ^Department of Life Sciences, Texas A&M 
University - Corpus Christi, Corpus Christi, TX 78412, USA, **Instituto Politecnico Nacional, Centra Interdisciplinary de 
Ciencias Marinas, Coleccibn Ictiologica, La Paz, Baja California Sur, 23096 Mexico 



We combine kinship estimates with traditional F-statistics to explain contemporary 
drivers of population genetic differentiation despite high gene flow. We investigate 
range-wide population genetic structure of the California spiny (or red rock) lobster 
(Panulirus interruptus) and find slight, but significant global population differentiation 
in mtDNA (* ST = 0.006, P = 0.001; D est chao = 0.025) and seven nuclear microsatellites 
(F ST = 0.004, P < 0.001; D esl C hao = 0.03), despite the species' 240- to 330-day pelagic lar- 
val duration. Significant population structure does not correlate with distance between 
sampling locations, and pairwise F ST between adjacent sites often exceeds that among 
geographically distant locations. This result would typically be interpreted as unex- 
plainable, chaotic genetic patchiness. However, kinship levels differ significantly 
among sites (pseudo-F 16 988 = 1.39, P = 0.001), and ten of 17 sample sites have signifi- 
cantly greater numbers of kin than expected by chance (P < 0.05). Moreover, a higher 
proportion of kin within sites strongly correlates with greater genetic differentiation 
among sites (D esl chao , R 2 = 0.66, P < 0.005). Sites with elevated mean kinship were 
geographically proximate to regions of high upwelling intensity (R 2 = 0.41, P = 0.0009). 
These results indicate that P. interruptus does not maintain a single homogenous pop- 
ulation, despite extreme dispersal potential. Instead, these lobsters appear to either 
have substantial localized recruitment or maintain planktonic larval cohesiveness 
whereby siblings more likely settle together than disperse across sites. More broadly, 
our results contribute to a growing number of studies showing that low F ST and high 
family structure across populations can coexist, illuminating the foundations of cryptic 
genetic patterns and the nature of marine dispersal. 
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Determining the temporal and spatial scales of dispersal 
and gene flow across species' ranges is essential for 



Introduction 



effective conservation and management. F-statistics 
(Wright 1943) and their analogues (e.g. Nei 1973; Weir 
& Cockerham 1984; Excoffier et al. 1992; Hedrick 1999) 
have been the workhorses in this regard for over 
65 years. However, as both the number and diversity of 
genetic markers have increased, so has the demand for 
analyses that can complement fixation indices and 
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extend our understanding of genetic data beyond the 
single marker, two-allele system pioneered by Wright 
(1943). Coalescent simulations (Kingman 1982) have 
emerged as the most informative techniques for distin- 
guishing between historical and contemporary drivers 
of population differentiation (Hey & Wakeley 1997; 
Tavare et al. 1997; Rosenberg & Nordborg 2002; Rozas 
et al. 2003; Drummond et al. 2005; Hickerson et al. 2006; 
Hey & Nielsen 2007; Eldon & Wakeley 2009). By incor- 
porating data from multiple nucleotide sequence-based 
markers, these equilibrium-independent analyses can 
isolate the effects of effective population size, demo- 
graphic history, migration, mutation and drift summa- 
rized by F ST (Hart & Marko 2010; Marko & Hart 2011, 
2012). However, for fragment length data such as that 
generated by microsatellite markers, a number of alter- 
native approaches have been advanced that can add 
insight into what is driving F ST patterns (reviewed in 
Hedgecock et al. 2007; Lowe & Allendorf 2010). 

One underutilized approach is the coupling of indi- 
rect metrics of gene flow (e.g. F-statistics, D esl chao ) with 
more direct measures such as kinship or parentage 
analyses (e.g. Loiselle et al. 1995; Selkoe et al. 2006; 
Buston et al. 2009; Christie et al. 2010; Palsboll et al. 
2010). Broadly speaking, kinship analyses provide an 
index of the relative relatedness of all genotyped indi- 
viduals in a data set, and parentage is a distinct case of 
kinship whereby the most likely parents of individual 
juveniles are identified (Vekemans & Hardy 2004; Jones 
& Arden 2003; reviewed in Blouin 2003; Jones et al. 
2010). Kinship coefficients (also known as coefficients of 
coancestry) are widely interpreted as the probability of 
identity by descent of the genes, but they are more 
properly defined as 'ratios of differences of probabilities 
of identity in state' (Hardy & Vekemans 2002, p. 23) 
from homologous genes sampled randomly from each 
pair of individuals (Hardy & Vekemans 2002; Rousset 
2002; Blouin 2003; Vekemans & Hardy 2004). 

By comparison, F-statistics and D est chao are often 
blind to the relatedness of individuals; different popula- 
tion samples with the same kinship structure can have 
very different levels of genetic differentiation among 
them and vice versa. By assessing how alleles are shared 
among individuals, kinship analyses can elucidate which 
locations have comparatively little ongoing genetic 
exchange in situations where low F$t values suggest 
high contemporary population connectivity. This clarifi- 
cation is important because such inferences can in fact 
be due to historically high migration rates, effective pop- 
ulation sizes or measurement error (Hart & Marko 2010; 
Marko & Hart 2011, 2012; Faurby & Barber 2012). 

Direct evidence of coancestry between individuals 
provides a particularly valuable complement to F-statis- 
tics when it is not possible to derive other independent 



estimates of demographic connectivity such as through 
the tagging and tracking of adults or larvae (e.g. Bell- 
quist et al. 2008; Meyer et al. 2009; Cartamil et al. 2011; 
Carson et al. 2011; Lopez-Duarte et al. 2012; reviewed in 
Lowe & Allendorf 2010). In marine systems, the major- 
ity of taxa have relatively sedentary adults, but a pela- 
gic larval stage that persists in the water column from a 
few minutes to over a year (Thorson 1950; Strathmann 
1987; McEdward 1995). These larvae are notoriously dif- 
ficult to track directly (Levin 2006), but the time that 
larvae spend in the open ocean has led to the intuitive 
expectation that the majority of marine species have 
high levels of gene flow (Hedgecock et al. 2007). How- 
ever, the preponderance of recent indirect genetic evi- 
dence, based mostly on F-statistics, indicates that there 
is generally a weak relationship between dispersal 
potential inferred from pelagic larval duration (PLD) 
and genetic structure (reviewed in Bradbury et al. 2008; 
Shanks 2009; Weersing & Toonen 2009; Riginos et al. 
2011; Selkoe & Toonen 2011). Furthermore, it is gener- 
ally overlooked that indirect gene flow via multigenera- 
tional stepping-stone dispersal at small scales can 
mimic direct gene flow across large scales (Puebla et al. 
2012). The relatively few studies that have directly mea- 
sured larval dispersal through larval tagging or parent- 
age analyses have bolstered the claim that many larvae 
have limited dispersal and often recruit back to their 
region of origin (Jones et al. 2005; Gerlach et al. 2007; 
Planes et al. 2009; Lopez-Duarte et al. 2012). Due to the 
nature of these direct (kinship/parentage) versus indi- 
rect (F-statistics) measures of population connectivity, it 
is possible that direct analyses may identify recruitment 
patterns that cannot be detected using traditional 
F-statistics (Waples & Gaggiotti 2006; Saenz-Agudelo 
et al. 2009; Palsboll et al. 2010). For example, Christie 
et al. (2010) found little genetic differentiation (max 
F ST = 0.0097) among populations of bicolor damselfish 
(Stegastes partitus) in Exuma Sound, Bahamas, but par- 
entage analysis identified high levels of self-recruitment 
at two of the eleven sampled locations. The direct iden- 
tification of parent-offspring pairs resulted in very dif- 
ferent management advice for this species than 
interpretation based on the F ST data alone. 

To date, most marine kinship studies have under- 
standably focused on parent-offspring identification in 
reef fishes with fairly short larval durations. Here, we 
show that kinship analyses can also be useful at the 
opposite end of the potential dispersal continuum: the 
California spiny, or red rock lobster, Panulirus interruptus 
(Randall 1840), spends at least the first 8 months of its 
life in the plankton, during which time it can presum- 
ably disperse across its entire geographic range. Species 
without barriers to dispersal are expected to exhibit no 
detectable neutral genetic population structure. Here, we 
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use the California spiny lobster as a model to test the 
intuitive assumption of genetic homogeneity in species 
with extended PLD. We demonstrate the utility of indi- 
vidual-based estimates of genetic exchange in the inter- 
pretation of connectivity based on F-statistics. 

Methods 

Study system 

The California spiny or red rock lobster, Panulirus inter- 
ruptus, exhibits high site fidelity during its adult phase 
(Withy-Allen 2010), but has a two-phase pelagic larval 
stage with a total PLD of 8-11 months (Johnson 1956, 
1960; Serfling & Ford 1975). The initial phyllosoma 
stage undergoes multiple moults to produce 11 verti- 
cally and geographically stratified stages in the pelagic 
environment (Johnson 1960; Pringle 1986). The 11th 
phyllosoma stage molts into the final puerulus stage, 
which settles into the lobster's preferred juvenile habi- 
tat. Panulirus interruptus can be found across a 1400-km 
Pacific coast range from Monterey Bay, CA (although 
very rare north of Point Concepcion) to Bahia Magda- 
lena, Mexico. Throughout its geographic distribution, 
P. interruptus plays an important ecological role in 
structuring both kelp forest and intertidal communities 
(Tegner & Levin 1983; Robles 1987; Lafferty 2004). 
Spiny lobsters are also a valuable commercial and recre- 
ational fisheries species in both Mexico and the USA 
with a combined value of over $39 million from the 
most recent estimates (Chavez & Gorostieta 2010; Porzio 
2012). 

Sample collection and DNA extraction 

We collected tissue samples from 1102 P. interruptus 
individuals across 17 sites located throughout the entire 
Pacific coastal range from Point Conception, CA, in the 
north to Bahia Magdalena (BMG), Baja California Sur, 
Mexico, in the south (Fig. 1, Table 1). Samples were col- 
lected nonlethally by removing a small piece of an 
antenna or a leg segment from each lobster. Lobsters 
were either captured by hand while Scuba diving or 
obtained from commercial fishermen. Tissue samples 
were preserved in 95% ethanol and stored at room tem- 
perature until extracted. DNA was isolated using a 
standard salting-out protocol (Sunnucks & Hales 1996), 
a rapid-boil technique (Valsecchi 1998) or DNeasy Ani- 
mal Tissue kits (Qiagen, Inc., Valencia, CA, USA). 

Mitochondrial DNA (mtDNA) 

We amplified a 494-bp fragment of cytochrome c oxi- 
dase subunit I gene (COI) using species-specific primers 



PintCOI-F (5'-GCTTGAGCTGGAATGGTAGG-3') and 
PintCOI-R (5'-CACTTCCTTCTTTGATCCC-3'), which 
were designed from GenBank sequence #AF339460 
(Ptacek et at. 2001) using primer3 (Rozen & Skaletsky 
2000). Polymerase chain reactions (PCRs) for each 
sample were performed in a 20-jj.I reaction volume con- 
taining 10 (.iL of 2x Biomix Red (Bioline, Taunton, MA, 
USA), 0.125 um each of forward and reverse primer, 
5-50 ng of genomic DNA and 0.75 x bovine serum 
albumin. PCR was carried out on a Bio-Rad Mycycler 
Thermal Cycler (Bio-Rad Laboratories Hercules, CA, 
USA), with an initial denaturation step of 95 °C for 
4 min, 35 cycles of denaturation (95 °C for 30 s), anneal- 
ing (56 °C for 30 s) and extension (72 °C for 30 s), fol- 
lowed by a final extension step of 72 °C for 10 min. 
PCR products were treated with 0.75 units of Exonucle- 
ase I and 0.5 units of Fast Alkaline Phosphatase (Exo- 
FAP; Thermo Fisher Scientific, Waltham, MA, USA) per 
7.5 |iL of PCR product and incubated at 37 °C for 
60 min, followed by deactivation at 85 °C for 10 min. 
Purified DNA fragments were sequenced in the 
forward direction with fluorescently labelled dideoxy 
terminators either on an ABI 3730XL capillary sequen- 
cer (Applied Biosystems, Foster City, CA, USA) by the 
Advanced Studies of Genomics, Proteomics and Bioin- 
formatics (ASGPB) Center at the University of Hawai'i 
at Manoa or an ABI 3130XL Genetic Analyzer (Applied 
Biosystems) at the Hawai'i Institute of Marine Biology 
EPSCoR Sequencing Facility. Unique sequences and 
sequences with ambiguous nucleotide calls were also 
sequenced in reverse to confirm sequence identity. 
Sequences were edited, aligned and trimmed to a uni- 
form size using sequencher 4.8b (GeneCodes Corpora- 
tion, Ann Arbor, MI, USA). The alignment did not 
contain any indels or frameshift mutations. 

We calculated nucleotide (tc) and haplotype diversity 
(h) for each sampling site as described in Nei (1987) 
using arlequin 3.5 (Excoffier et al. 2005). To visualize 
relationships between individual sequences, we con- 
structed a median-joining network (Bandelt et al. 1999) 
using network 4.6.0.0 (http://www.Fluxus-engineering. 
com/network terms.htm). We investigated population 
structure using an analysis of molecular variance 
(amova) as implemented in arlequin. We used an ana- 
logue of Wright's F S t (^st); which incorporates a model 
of sequence evolution, for both our complete data set and 
pairwise population comparisons (Weir & Cockerham 
1984; Excoffier et al. 1992). Using jModelTest2 (Guindon 
& Gascuel 2003; Darriba et al. 2012), we determined that 
the Tamura & Nei (1993) with a Ti/Tv ratio of 11.2 and 
gamma parameter of 2.1 was the most appropriate 
model of sequence evolution implemented in arlequin. 
Global $st and each pairwise population $ ST were 
tested for significance with 100,000 permutations. Due 
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Fig. 1 Map of lobster specimen collection 
locations in the Southern California Bight 
and Baja California, Mexico. CARP, Car- 
pinteria, CA; SMI, San Miguel Island, 
CA; SCI, Santa Cruz Island, CA; MLBU, 
Malibu, CA; SCAT, Santa Catalina Island, 
CA; SNI, San Nicholas Island, CA; 
CLEM, San Clemente Island, CA; CRDO, 
Islas Coronados, Mx; PTN, Puerto 
Nuevo, Mx; PBDA, Punta Banda, Mx; 
PBI, Punta Baja, Mx; IGP, Isla Guada- 
lupe, Mx; ODL, Laguna Ojo de Liebre, 
Mx; PEU, Punta Eugenia, Mx; BTG, Bahi- 
a Tortugas, Mx; ABRE, Punta Abreojos, 
Mx; BMG, Bahia Magdalena, Mx. Photo 
credit: Patrick W. Robinson. 



to the high heterozygosity in cytochrome c oxidase sub- 
unit I gene (COI) in P. interruptus, we also calculated 
D est chao as an absolute measure of differentiation 
between sites. The magnitude of F ST is inversely pro- 
portional to heterozygosity (Hedrick 2005; Meirmans 
2006; Jost 2008), while D est chao is less susceptible to 
biases caused by genetic diversity (Bird et al. 2011). For 
mtDNA, D est 

chao (Jost 2008) was calculated with spade 
(Chao & Shen 2009). Mantel tests using linearized F ST 
[Fst/ (1 — Fst)] and the natural log of Euclidean dis- 
tance (Rousset 1997) were conducted in genodfve 2.0b20 
(Meirmans & Van Tienderen 2004) to test for isolation 
by distance (IBD) among all sample sites, as well as 
subsets of sites across the range (California sites, Mex- 
ico sites, all island sites, all mainland / continental sites). 

Microsatellite DNA 

Seven previously developed microsatellites (A5, A102r, 
A110, Pin29L, Pinl89, PinlO and Pin244) (Ben-Horin 
et al. 2009) were amplified by PCR in a 10-|iL reaction 
volume containing 1 x GoTaq Reaction Buffer (with 
1.5 mM MgCl 2 pH 8.5; Promega Corp., Madison, WI, 
USA), 2 um total dNTPs, 0.1 U GoTaq polymerase 



(Promega Corp.), 6 |im each of forward and reverse pri- 
mer and 20 ng of genomic DNA. PCR was carried out 
using a Bio-Rad DNA Engine Dyad Thermal Cycler 
with the following conditions: an initial denaturation 
step of 95 °C for 3 min, 35 cycles of denaturation (94 °C 
for 40 s), annealing at primer-specific annealing temper- 
atures (for 40 s; see Table 1 in Ben-Horin et al. 2009) 
and extension (72 °C for 40 s) and a final extension step 
of 72 °C for 30 min. Forward primers were fluorescent- 
ly labelled with WellRed D2, D3 or D4 dye (Beckman 
Coulter Inc., Fullerton, CA, USA; see Table 1 in Ben- 
Horin et al. 2009 for primer labels). Microsatellite PCR 
products were sized on a Beckman Coulter CEQ 8000 
capillary sequencer with a 400-bp size standard (Beck- 
man Coulter Inc.). Alleles were scored using a CEQ 
8000 genetic analysis system (Beckman Coulter Inc.). 

Microsatellite quality control [Hardy-Weinberg equi- 
librium (HWE), linkage equilibrium, scoring errors] fol- 
lowed Selkoe & Toonen (2006) as detailed in Ben-Horin 
et al. (2009), but expanded across all sample locations 
and loci. Additionally, null allele frequency was re-cal- 
culated with ML-Relate (Kalinowski et al. 2006) to enable 
significance testing. We calculated observed (H Q ) and 
expected (H e ) heterozygosities at each location, as well 
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Table 1 Summary statistics for Panulirus interruptus listed from northernmost to southernmost collection sites: total number of indi- 
viduals sequenced for seven microsatellites (N, nDNA) and mtDNA cytochrome c oxidase subunit I (N, mtDNA) 
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Islas Coronados (CRDO) 
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Isla Guadalupe (IGP) 
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Laguna Ojo de Liebre (ODD 
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Punta Eugenia (PEU) 
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Bahia Tortugas (BTG) 
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Punta Abreojos (ABRE) 
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For mtDNA: haplotype diversity (h), effective number of haplotypes (h eSf ) and nucleotide diversity (ti). For nDNA microsatellites: rar- 
efied allelic richness (AR), effective number of alleles (A e ff), observed (H„) and expected (H e ) heterozygosity. 



as allelic richness (A), rarefied to 25 individuals, using 
the Excel Microsatellite Toolkit 3.1.1 (Park 2001) and 
fstat (Goudet 2001), respectively. We also calculated the 
effective number of alleles at each locus in GenoDive. 
Overall genetic structure was analysed using an amova 
framework (Excoffier et al. 1992) implemented in geno- 
dive assuming the infinite allele model, with Fgx equiva- 
lent to Weir & Cockerham's 9 (1984). Consistent with 
mtDNA (COI), microsatellite allelic diversity was like- 
wise high, so we calculated D est chao (Jost 2008) in addi- 
tion to F ST . GenoDive was used to calculate pairwise 
F ST and D est chao for all sampling location pairs, and 
F ST was tested for significance using 100 000 permuta- 
tions. Patterns of IBD were investigated for microsatel- 
lite data in GenoDive as described above for mtDNA. 
The program geste 2.0 (Foil & Gaggiotti 2006) was used 
to calculate local F ST , a site-specific metric of allelic dif- 
ferentiation that accounts for the nonindependence 
inherent in multiple comparisons (Balding & Nichols 
1995; Hudson 1998). 

Kinship 

In order to further examine the factors driving the 
observed F sx and D est chao patterns, we calculated kin- 
ship coefficients (Loiselle et al. 1995) for each pair of 
individuals in genodive. These coefficients are based on 
the probability of identity of two alleles for each pair of 
homologous genes compared between each pair of 



individuals. Kinship was estimated with respect to the 
allele frequencies for the full data set, so these coeffi- 
cients provide an index of relative relatedness between 
each pair of individuals. In order to determine whether 
individuals collected at the same location were more 
closely related to each other than individuals collected 
at different locations, we conducted an amova on the 
kinship coefficients. This approach compared the varia- 
tion in within-population kinship coefficients with the 
variation in among-population kinship coefficients 
using the permanova+ 1.0.2 software add-on as imple- 
mented in primer6 (Clarke & Warwick 2001), following 
Stat et al. (2011) and Padilla-Gamiho et al. (2012). Specif- 
ically, the kinship covariance matrix created in genodive 
was loaded into primer6 as a Correlation Resemblance 
Matrix. A simple one-way analysis of variance (termed 
permanova in primer6) was conducted with 10 000 unre- 
stricted permutations of the raw data and type III sums 
of squared differences. Significance is determined by 
evaluating a pseudo-F value (Clarke & Warwick 2001) 
based on the F-distribution, which is not to be confused 
with Wright's F-statistics. 

To investigate site-specific patterns in kinship, we 
counted the number of closely related individuals 
within each site where the kinship coefficient was 
greater than or equal to the equivalent of that expected 
for quarter-siblings (0.047). Following Selkoe et al.'s 
(2006) and Buston et al.'s (2009) analysis of relatedness 
approach, we binned our counts of closely related 
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lobsters according to specific levels of kinship. We used 
the Loiselle et al.'s (1995) coancestry coefficients (full- 
sib = 0.25, half-sib = 0.125) to generate the following 
bins: 'nearly identical', 0.57 > k > 0.375; 'full-sib', 
0.375 > k > 0.1875; 'half-sib', 0.1875 > k > 0.09375; and 
'quarter-sib', 0.09375 > k > 0.047. These bounds repre- 
sent the midpoints between the coancestry coefficients 
in Loiselle et al. (1995), and we use quarter-sib as a 
short-hand to represent half of the level of coancestry 
as half-sibs. The nearly identical bin represents compar- 
isons above full-sibs and is based on our kinship coeffi- 
cient distribution for comparisons of individuals to 
themselves. We tested multiple bin sizes and divisions 
and found our results to be robust to these changes. To 
test for an overabundance of closely related lobsters 
within sites, we implemented a permutation test (10 000 
replicates) where the lobsters were randomly assigned 
to sites and the observed number of closely related 
individuals was compared to the null distribution for 
each site. We calculated the maximum-likelihood esti- 
mate of relatedness (r) in ML-Relate, following the scale 
for the index of relatedness (full-sibs, r = 0.5; half- 
sibs = 0.25) to be able to compare our kinship results 
with a relatedness index. We also tested the relationship 
of the mean pairwise F ST , mean pairwise D e st_chao and 
local F ST for each site with the proportion of closely 
related lobsters at that site (summed across all four sib- 
ship categories for the Loiselle et al. (1995) metric and 
across half-sibs and full-sibs in ML-Relate). 

Upivelling 

We identified hotspots of genetic differentiation by cal- 
culating mean D est chao for each sampling location. We 
interpolated these values beyond the specific sites we 
sampled using an inverse weighted distance (IWD) 
algorithm in the Spatial Analyst extension in arcgis 10. 
To test the hypothesis that upwelling is a potential dri- 
ver of both increased kinship and, in turn, site-specific 
genetic structure for P. interruptus in this region, we 
overlaid on this map known areas of consistent upwell- 
ing in Baja California, Mexico (as identified by Zaytsev 
et al. (2003)). We tested the relationship between the 
mean kinship at a site and that site's closest distance to 
an area of persistent upwelling. For sites within the 
Southern California Bight (SCB), where there are no 
areas of persistent upwelling, either the distance to 
Point Conception or to the edge of the northernmost 
upwelling area in Baja California, Mexico, was used, 
whichever was shorter. The southernmost sites in the 
range were measured to an upwelling area just off of 
BMG, which was identified by Zaytsev et al. (2003), but 
is not included in our figure. Negative distances repre- 
sent sites that are located within upwelling regions and 



are measured from their location to the nearest edge of 
the upwelling zone. Both kinship and upwelling regres- 
sion analyses were performed in spss 17.0. 

Results 

Mitochondrial DNA (mtDNA) 

We sequenced COI for 931 individuals across 17 sites, 
which yielded 238 haplotypes. Haplotype diversity (h) 
was high, ranging from 0.88 to 0.95 (8.3-20 effective 
haplotypes), with a mean of 0.92 (12.5 effective haplo- 
types). In contrast, nucleotide diversity (rc) was rela- 
tively low, ranging from 0.005 to 0.018, with a mean of 
0.010. The number of individuals sequenced (N), haplo- 
type diversity (h), effective number of haplotypes (h e(f ) 
and nucleotide diversity (n) for each site are listed in 
Table 1. 

The median-joining network (Fig. 2) reveals two 
dominant haplotypes differing by one-nucleotide substi- 
tute and present at all sample sites. The most numerous 
haplotype was found in 235 individuals (25% of indi- 
viduals sequenced), and the second most numerous 




Fig. 2 Median-joining network for Panulirus interruptus 
mtDNA, constructed using 454 base pairs of cytochrome c oxi- 
dase subunit I (COI) from each of 931 individuals in the pro- 
gram network 4.6.0.0. Each circle is a unique haplotype 
proportional in size to the number of individuals with that 
haplotype. The two largest circles represent 235 and 95 indi- 
viduals. The smallest circle represents two individuals: there 
are 131 singletons in the data set, but these have been omitted 
for ease of visualization. A full, unedited network is included 
in the supplementary material (Fig. S3, Supporting informa- 
tion). Colours correspond to one of 17 locations where the indi- 
vidual haplotypes were found (see key, Fig. 1, Table 1). Lines 
connecting haplotypes represent a single base-pair difference 
between haplotypes, with crossing lines representing addi- 
tional differences. 
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haplotype was found in 92 individuals (10%). The hap- 
lotype network was characterized by a starburst pat- 
tern, with the majority of remaining haplotypes 
differentiated by one to two base pairs from the domi- 
nant haplotypes. Eighty haplotypes were represented in 
only two individuals and there were 131 singletons 
across all 17 sites; the removal of singletons did not 
impact the overall structure of the haplotype network 
(Fig. 2), so they were omitted for ease of visualization. 
A full, unedited network is included in Fig. S3 (Sup- 
porting information). 

Global $st (0.006) was low, but significant 
(P = 0.001). Statistically significant pairwise <& ST values 
were 0.01-0.04, and in general, D est chao values were an 
order of magnitude higher than $ ST , both for global 
D e st_chao (0.025) and for each of the pairwise compari- 
sons, ranging up to 0.300 (Table 2). Among pairwise 
$ sx comparisons, 26 of 136 (19%, Table 2) were signifi- 
cant at P < 0.05; however, after correcting for the false 
discovery rate (Benjamini et al. 2006), no pairwise com- 
parisons were significant (P < 0.00035, Table 2). There 
was no pattern of IBD (Wright 1943) at any scale in the 
data set. This was true whether we examined all sam- 
pling locations together, or for any of the specific subre- 
gions: all Mexico sampling locations, California 
locations (SCB), only islands and only continental sites 
(Table SI, Supporting information). 

Microsatellite DNA 

We scored 989 individuals across 17 sites for seven 
nuclear microsatellite loci, with five to 104 alleles per 
locus, which translated to between two and 23 effective 
alleles per locus (Table S2, Supporting information). 
The genotyping error rate, determined by re-genotyping 
74 individuals at each locus, ranged from 0.0% to 4.1%, 
with an average of 2.3% overall (Table S2, Supporting 
information). We identified two pairs of individuals 
with identical genotypes and removed one individual 
of each pair from the data set in order to eliminate the 
possibility that the same individual was sampled two 
times. In both cases, identical genotypes were identified 
within a sampling site. The expected chances of observ- 
ing true identical twins in this study ranged from 
1 x 10~ 42 to 4.39 x 10 -26 , while our observed rate was 
substantially greater, 2 in 989 specimens. 

Rarefied allelic richness was similar among sites and 
ranged from 16.19 to 18.51, while the effective number 
of alleles ranged from 12.05 to 15.98 (Table 1). Expected 
heterozygosity (ff e ) was between 0.86 and 0.90, while 
observed heterozygosity (H D ) exhibited a slightly wider 
range from 0.76 to 0.87 (Table 1). Tests for linkage dis- 
equilibrium (LD) were significant in 20 of 349 compari- 
sons (-6%) after correcting for multiple tests, and there 



were no locus-specific patterns. The three sample sites 
with the highest percentage of LD comparisons corre- 
spond to the three sites with the highest levels of 
kinship, as reported below, suggesting kinship may be 
high enough at these sites to produce a signal of LD. 
In general, however, LD is a weak test of family struc- 
ture, and we found no site-specific patterns across the 
rest of the sites in the study. There were significant 
deviations from HWE in 46 of 119 (-39%) comparisons 
after correcting for multiple comparisons, but again no 
sample-specific patterns were observed, micro-checker 
2.2.3 (Van Oosterhout et al. 2004) found no evidence of 
scoring errors due to large allele dropout or stutter; 
however, six of seven markers showed patterns consis- 
tent with null alleles, which are the likely cause of the 
deviations from HWE (see Ben-Horin et al. 2009). 
Although the frequencies of these null alleles are low 
(1.20-6.63% across loci, Table S2, Supporting informa- 
tion), we wanted to be sure they would not affect our 
results. To test for the impact of null alleles on our 
results, we used FreeNA (Chapuis & Estoup 2007) to 
generate alleles for the data set where nulls were 
expected, re-analysed the data, and our subsequent 
results and conclusions remained the same. Therefore, 
we present only analyses with the full original data set. 

Significant partitioning of the samples among loca- 
tions was detected in the microsatellite loci with both a 
global fixation test (F ST = 0.004, P < 0.001) and a global 
genetic differentiation test (D estChao = 0.03, P < 0.0005). 
Global partitioning was also detected in analyses of 
individual loci, with a significant global F ST (0.002- 
0.011, Table S2, Supporting information) and global 
D es t_chao (0.030-0.114) at each individual locus, except 
A110. 

We ran pairwise F sr and D est chao comparisons for 
each marker individually, as well as jackknifing across 
markers, and found the results were fairly consistent 
throughout these comparisons. Significant pairwise Pst 
comparisons among sampling sites using the microsat- 
ellite loci were fairly low, ranging from 0.002 to 0.015, 
but 71 of 136 comparisons (52%, Table 2) remained sig- 
nificant after correcting for the false discovery rate 
(Benjamini et al. 2006). As in our mtDNA results, pair- 
wise D est _chao comparisons generally were an order of 
magnitude higher than each respective pairwise F ST 
comparison and ranged from —0.021 to 0.128 (Table 2). 
The higher magnitude of D est chao compared F ST to 
matches the expectation when heterozygosity is -0.9, as 
it is in this study (Bird et al. 2011). 

The sites with the highest mean pairwise P ST and 
D est chao values were in the northern and central Baja 
California, Mexico region [Puerto Nuevo (PTN), Punta 
Banda (PBDA), Punta Baja and Bahia Tortugas (BTG)]. 
These sites also stood out by having the highest 
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proportion of significant differences when compared to 
the other sites. Three of these four sites also had the 
highest local F ST values (PBDA, BTG and PTN). There 
was no pattern of IBD in the data set, whether we 
examined all sampling locations together, or for any of 
the specific subregions, as described for the mtDNA 
results (Table SI, Supporting information). 

Kinship 

Following the exclusion of lobster specimens with identi- 
cal genotypes, as well as self-comparisons, kinship coef- 
ficients ranged from —0.155 to 0.57 (Fig. 3). The overall 
mean kinship, which is expected to be zero, was 
-0.000025 ± 0.00008SE. There was a disparity, however, 
between the within-site mean kinship (0.003 ± 0.0004) 
and the among-site mean kinship (-0.0002 ± 0.00009) 
(Fig. 3). Mean within-site maximum-likelihood estimates 
of relatedness (r, from ML-Relate) ranged from 0.034 to 
0.057 (Table S3, Supporting information) and were sig- 
nificantly correlated with mean kinship values (ft:) at each 
site (Table S3, Supporting information, Pearson's 
r = 0.92, P < 0.0005). Kinship coefficients were signifi- 
cantly greater for within-site than for among-site com- 
parisons (pseudo-Fi6,988 = 1-39, P = 0.001). In total, 10 of 
17 sites had significantly greater numbers of closely 
related pairs of individuals than expected by chance 
(P < 0.05) in at least one of four kinship categories: five 
sites had an excess of individuals in the 'nearly identical' 
category (0.57 > k > 0.375), five sites had an excess in the 
'full-sib' group (0.375 > k > 0.1875), five sites had an 
excess in the 'half-sib' category (0.1875 > k > 0.09375), 
and four sites had an excess in the 'quarter-sib' bin 
(0.09375 > k > 0.047) (Fig. 4, Table S3, Supporting infor- 
mation). The proportion of kin in each site was signifi- 
cantly related to mean pairwise F ST (R 2 = 0.669; 
F U6 = 30.333, P < 0.0005) and D est C h ao (R 2 = 0.658; 
F U6 = 28.885, P < 0.0005; Fig. 5) for each site as well as 
to local F ST (R 2 = 0.243; F 1/16 = 4.825, P = 0.044; Fig. 5). 
These findings are consistent across kinship classes. 
When we removed the lowest level of kinship (examin- 
ing only kinship levels equivalent to 'nearly identical', 
'full-sib' and 'half-sib'), the relationship of both F ST 
(R 2 = 0.672; F W6 = 30.741, P < 0.0005) and D est C hao 
(R 2 = 0.658; F 146 = 28.887, P < 0.0005) with the propor- 
tion of kin stayed approximately the same, while the 
relationship of local Fst with the proportion of kin 
strengthened (R 2 = 0.303; F 1/16 = 6.526, P = 0.022). The 
relationships of mean pairwise F ST , D est chao and local 
F sx with the proportion of related individuals (r, full- 
and half-sibs) from ML-Relate were all significant and 
stronger than the relationships between these summary 
statistics and proportion of kin (Fig. SI, Supporting 
information). 



Upwelling 

Regions containing persistent, strong upwelling 
regimes, both across seasons within a year and across 
all years from 1996 to 2002, were traced from fig. 14 in 
Zaytsev et al. (2003) and are depicted by dotted lines in 
Fig. 6a. Three of the four upwelling regions overlap the 
geographic range of our study and contain four of our 
sampling locations: Islas Coronados, PTN, PBDA and 
BTG (Fig. 6a). Three of these four sites (PTN, PBDA 
and BTG) contain one of the highest four values for 
mean D est chao , mean local Fst and mean kinship. The 
IWD function applied to this data extrapolates these 
site-specific metrics across the entire study area, includ- 
ing the unintended extrapolation of D est chao values to 
ocean areas where adult lobsters do not live, but 
phyllosoma may be present. Along the Baja California 
coastal areas containing adult lobsters, all three of the 
areas with the highest genetic differentiation (D est chao) 
overlap the regions of strong upwelling intensity 
(Fig. 6a). The results for mean kinship and local Fst 
were similar, although not shown. Across the geo- 
graphic extent of our study, we see a significant rela- 
tionship (R 2 = 0.407; F 1/16 = 10.296, P = 0.006) between 
mean kinship at each site and the closest distance 
between each site and the edge of an upwelling zone 
[(mean kinship + 1) = —0.005 In (distance to upwell- 
ing + 100 km) + 1.028, Fig. 6b). This relationship holds 
for mean relatedness (r) versus distance from upwelling 
as well (Fig. S2, Supporting information). 

Discussion 

There is accumulating evidence from multiple 
approaches that larvae rarely reach their full dispersal 
potential, resulting in a paradigm shift away from the 
perception that most marine populations are genetically 
homogenous across broad geographic scales (Jones et al. 
1999; Swearer et al. 1999, 2002; Mora & Sale 2002; Gran- 
tham et al. 2003; Taylor & Hellberg 2003; Marko 2004; 
Cowen et al. 2006; Becker et al. 2007; Lopez-Duarte et al. 
2012). However, this evidence has come exclusively 
from species with short to modest larval periods 
(1-60 days). Even among species with modest PLD, 
there are just a few striking examples of broad popula- 
tions with no genetic substructure across their full 
range. For example, reef fishes Myripristis jacobis in the 
Atlantic Ocean (Bowen et al. 2006), Lutjanus kasmira in 
the Central Pacific and Eastern Indian Oceans (Gaither 
et al. 2010) and Acanthums nigrofuscus in the Pacific 
(Eble et al. 2011) all exhibit genetic homogeneity across 
thousands of kilometres (up to 12 000 km). Dispersal 
potential is assumed to be great in species with very 
long PLDs (>120 days), and population genetic surveys 
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Fig. 3 (a) Line graph depicting the total 
number of kinship (k) pairwise compari- 
sons for each 0.01 bin of kinship (solid 
line, total N = 457,652) and the number 
of within-site kinship comparisons 
(dashed line, total N = 30,914). (b) Inset 
of (a), depicting the total number of kin- 
ship comparisons (solid line) and the 
number of within-site comparisons 
(dashed line) from 0.25 to 0.57 on a sepa- 
rate j/-axis that ranges from 0 to 140 com- 
parisons. The majority of these kinship 
comparisons are between individuals 
sampled at the same location, (c) Distri- 
bution of kinship coefficients divided 
into 0.01 bins and coloured by the pro- 
portion of within-site (dark grey) versus 
among-site (light grey) comparisons 
within each 0.01 division. White bars rep- 
resent levels of kinship that were not 
found in the data set. Bars on the *-axis 
represent the divisions between unre- 
lated and related individuals and 
between each of the four kinship catego- 
ries we analysed: 'quarter-sib', 0.047 < 
k < 0.09375; 'half-sib', 0.09375 < k < 
0.1875; 'full-sib', 0.1875 < k < 0.375; and 
'nearly identical', 0.375 < k < 0.57. 



Nearly Identical 



Kinship Coefficient 



of such species to date have revealed little population 
structuring across broad geographic scales (Ovenden 
et al. 1992; Silberman et al. 1994; Thompson et al. 1996; 
Tolley et al. 2005; Garcia-Rodriguez & Perez-Enriquez 
2006; Inoue et al. 2007; Home et al. 2008; Reece et al. 
2011). When genetic discontinuities have been observed 
in species with long PLD, they have invariably corre- 
sponded with known biogeographic barriers, or oceano- 
graphic transitions (Palero et al. 2008; Babbucci et al. 
2010; Chow et al. 2011). 

In contrast to the intuitive expectation that Panulirus 
interruptus, with a minimum PLD of 240 days, would 
be genetically homogenous across its entire 1400-km 
range along the west coast of North America, we found 
slight, but significant genetic structuring among several 
sampling locations throughout Mexico and Southern 
California (Table 2). This finding contrasts with previ- 
ous work that did not detect population structure in 



P. interruptus throughout Baja California, Mexico, using 
mtDNA RFLPs (Garci'a-Rodriguez & Perez-Enriquez 
2006). Notably, lobsters do not exhibit a genetic break 
across Punta Eugenia, a faunal boundary for rocky 
intertidal species (Valentine 1966; Blanchette et al. 2008; 
Gaines et al. 2009) and a phylogenetic break for a num- 
ber of coastal fishes (Bernardi et al. 2003). Nor does the 
overall pattern of genetic differentiation in P. interruptus 
correspond to the Northern, Central and Southern 
regional population subdivision within Baja predicted 
by Perez-Enriquez et al. (2001). Rather, genetically dif- 
ferentiated sites are nested within a greater area of 
undifferentiated sites (Fig. 6a). Specifically, some sites 
exhibit no genetic differentiation across the 1400-km 
species range, whereas other sites are differentiated 
from the majority of sampled sites, and there is no sig- 
nal of IBD in either the mtDNA or the nuclear microsat- 
ellite markers across multiple spatial scales (Table SI, 
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Fig. 4 Histogram of the proportion of kinship observed for each site that is in excess of expected levels due to chance in four kinship 
categories: 'quarter-sib', 0.047 < k < 0.09375; 'half-sib', 0.09375 < k < 0.1875; 'full-sib', 0.1875 <k< 0.375; and 'nearly identical', 
0.375 < k < 0.57. Asterisks indicate significant (P < 0.05) differences between the observed and expected proportion of kinship com- 
parisons within that category at that site. Expected kinship levels were constructed using 10 000 permutations of all kinship values 
across all sites without replacement to generate the distribution of kinship values that should occur if individuals were randomly dis- 
tributed among sites. The nearly identical bin represents comparisons above full-sibs and is based on our kinship coefficient distribu- 
tion for comparisons of individuals to themselves. Counts of the actual number of pairwise comparisons at each site that fell within 
each kinship category are listed in Table S3. 
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Fig. 5 Linear regression of mean pairwise D eal chao (circles, 
solid line) and local F sx (squares, dashed line) for seven micro- 
satellite loci. Both metrics were regressed on the proportion of 
within-site kinship comparisons at each site that are signifi- 
cantly (P < 0.05) greater than k = 0.047 (the lower boundary of 
the four designated kinship categories). 

Supporting information). Similar patterns of genetic dif- 
ferentiation among proximate sites have been shown in 
species with shorter larval developmental periods 
(reviews by Larson & Julian 1999; Hedgecock et al. 
2007; Riginos et al. 2011; Toonen & Grosberg 2011), but 
have not been reported for a species with such a long 



PLD as this one. Additionally, P. interruptus has equiva- 
lent or greater levels of genetic substructure than other 
species occurring in this same region (Selkoe et al. 

2010) , despite a PLD that is an order of magnitude 
higher. 

The term 'chaotic genetic patchiness' was coined 
(Johnson & Black 1982, 1984) to describe ephemeral, 
finely spatio-temporal patterns of genetic structure gen- 
erated by variation in the larval pool, recruitment and 
natural selection, which are counteracted in the long- 
term by dispersal and gene flow (Toonen & Grosberg 

2011) . Much of the difficulty in interpreting these unex- 
pected patterns in genetic differentiation is due to the 
nature of Fg T as a summary statistic. Significant struc- 
ture among populations may be a result of differences 
in effective population size (and corresponding genetic 
drift), demographic or colonization history, migration 
or some combination of these factors, especially for 
populations that may not have reached migration-drift 
equilibrium. Direct interpretation of summary statistics 
(F ST , DesLChao) in the context of gene flow can be prob- 
lematic (reviewed by Lowe & Allendorf 2010; Hart & 
Marko 2010; Marko & Hart 2011; Bird et al. 2011; Karl 
et al. 2012), especially in species such as P. interruptus, 
with highly fecund individuals and a potential for 
reproductive skew (Eldon & Wakeley 2009). For the 
many marine species with high fecundity and a type III 



© 2013 John Wiley & Sons Ltd 



SITE- SPECIFIC KINSHIP IN THE RED ROCK LOBSTER 3487 




Fig. 6 (a) Map of lobster specimen collection locations overlaid 
with an inverse weighted distance interpolation of mean 
D es t_chao values at each site created using the Spatial Analyst 
extension in arcgis 10. Red coloration represents areas of great- 
est genetic differentiation, green represents areas of little 
genetic differentiation, and yellow and orange represent the 
graded values between these extremes. Mean kinship values at 
each site showed the same pattern as D est chRO , with the high- 
est kinship in red and lowest in green. Dashed lines circle 
areas of consistently high upwelling intensity (adapted from 
Zaytsev et al. 2003). There are also areas of high upwelling 
intensity at Point Conception, and one south of Bahia Magda- 
lena that are not captured by this map. (b) Log-linear regres- 
sion of mean kinship at each site on the distance (km) to the 
nearest edge of an area of high upwelling intensity [from (a)]. 

survivorship curve, an independent test can help deter- 
mine whether gene flow is the primary driver resulting 
in the observed population structure (Hart & Marko 
2010; Lowe & Allendorf 2010). 

Here, kinship (Loiselle et al. 1995) enriches our under- 
standing of the drivers underlying significant 
differences in F ST among sites. The pattern of chaotic 
genetic patchiness in P. interruptus evident in the F ST 
and D est _ chao analyses seems to be primarily a result of 
the nonrandom occurrence of closely related lobsters 
within sample sites. Across all sites, lobsters were more 
closely related within sites than between sites, and at 
the majority of sites, we found significantly greater than 
expected levels of kinship between adult lobsters 
(Fig. 4). Moreover, the proportion of kin found at each 
site accounts for the majority of the variation in the 
sites' genetic differentiation: the most greatly 



differentiated sampling sites have the highest propor- 
tion of kin (Fig. 5). 

One potential scenario that could generate high pro- 
portions of kin within sites is recruitment pulses of 
related individuals. The simplest explanation for this 
phenomenon is that larvae released together stay 
together throughout dispersal and recruitment (kin 
aggregation). However, this pattern would also result 
from extreme differential reproductive success among 
individuals, so that a recruiting cohort is entirely made 
up of offspring from only a few individuals (sweep- 
stakes recruitment, Hedgecock 1986, 1994a,b). Sweep- 
stakes recruitment could also generate the star-shaped 
pattern of our mtDNA haplotype network (Fig. 2, Fig. 
S3, Supporting information), although this pattern could 
also be indicative of a recent population expansion. Pre- 
vious kinship analyses have detected high levels of 
relatedness within cohorts of larval recruits in both 
fishes (Planes et al. 2002; Pujolar et al. 2006; Selkoe et al. 
2006; Buston et al. 2009; Bernardi et al. 2012) and inver- 
tebrates (Veliz et al. 2006), which supports both the 
hypothesis of kin aggregation throughout development 
and/or the hypothesis of sweepstakes reproduction. 
Unfortunately, we could not directly test these alterna- 
tives in P. interruptus because we did not have samples 
of new recruits. Nevertheless, given the size selectivity 
of lobster traps, and the intense fishing pressure for lob- 
ster depressing the age range (Iacchei et al. 2005; Kay & 
Wilson 2012), it is possible that our samples are largely 
made up of single year classes consisting of closely 
related individuals recruiting together by one of these 
aforementioned mechanisms. To our knowledge, this 
study is the first documented case of kin aggregation in 
the adult population of a marine species with plank- 
tonic larvae, although kin aggregation in recruits has 
been reported. Previous studies that have looked at 
only kin relationships among adults, rather than among 
cohorts of recruits, have found no evidence of kin 
aggregation in marine species (Avise & Shapiro 1986; 
Kolm et al. 2005; Buston et al. 2007; Palm et al. 2008; 
Andrews et al. 2010; Berry et al. 2012). Consequently, 
kin aggregation is generally assumed to be a transient 
phenomenon limited to newly settled recruits, with lit- 
tle detectable signal in adult populations due to multi- 
ple source populations of recruits, changes in 
reproductive success and differential juvenile mortality 
(Kordos & Burton 1993; Moberg & Burton 2000; Flowers 
et al. 2002; Planes et al. 2002; Selkoe et al. 2006; Buston 
et al. 2009). 

High levels of within-site kinship could also be driven 
by a temporally stable pattern of self-recruitment, either 
through larval retention or through larval dispersal with 
subsequent recruitment back to the natal site. The pros- 
pect that larvae stay in the plankton for 240-330 days 
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and return to settle near their site of release seems unli- 
kely at first. However, the site-specific kinship patterns 
in our data match theoretical predictions for a species 
that has evolved a long PLD to avoid predation during 
the larval phase rather than to facilitate broad dispersal 
of larvae (Strathmann et al. 2002). The extended PLD 
may enable phyllosoma to disperse far offshore, into a 
pelagic environment that is favorable for the survival of 
unprotected larval-stage individuals (Strathmann et al. 
2002). Late-stage lobster larvae (pueruli) are fast swim- 
mers (Serfling & Ford 1975) and may utilize strong 
upwelling regimes to return and settle near their natal 
site after dispersing offshore. If this behaviour is selec- 
tively advantageous, we would expect to observe 
enhanced local recruitment regardless of PLD. Further- 
more, local recruitment should be more pronounced at 
sites with strong, persistent upwelling (Fig. 6b). 

The four sites with the highest proportion of kin in 
Baja California, Mexico, are located in areas of persis- 
tent upwelling (Zaytsev et al. 2003; Fig. 6a). Although 
upwelling was initially proposed as a mechanism for 
advecting recruits away from coastal areas (Roughgar- 
den et al. 1988), subsequent studies have questioned 
that prediction (e.g. Shanks & Brink 2005; Morgan et al. 
2009; Morgan & Fisher 2010; reviewed in Shanks & 
Shearman 2009). Genetic studies have proven equivocal 
in this regard: some have found reduced genetic sub- 
structure between populations in years with greater 
cumulative upwelling (e.g. Flowers et al. 2002; Barshis 
et al. 2011), and others have found upwelling dynamics 
to be insufficient to explain temporal and spatial genetic 
patterns (e.g. Toonen & Grosberg 2011). Passive larval 
dispersal models indicate that the effects of upwelling 
may depend on whether larvae stay near the surface 
(upwelling advects larvae offshore) or undergo regular 
migrations to depth (upwelling delivers larvae to the 
coast) (Byers & Pringle 2006; Marta-Almeida et al. 2006). 

Spiny lobsters are known to have relatively large lar- 
vae capable of dynamic movement. Both the phylloso- 
ma and puerulus stages show evidence of active 
movement in the pelagos, with phyllosoma exhibiting 
diel vertical migrations as well as horizontal move- 
ments (Kittaka 1994; Chiswell & Booth 1999; Phillips 
et al. 2006; Butler et al. 2011), and pueruli demonstrating 
rapid swimming, navigation towards the coast and hab- 
itat settlement preferences (Serfling & Ford 1975; Jeffs 
et al. 2005; Phillips et al. 2006). In P. interruptus specifi- 
cally, Pringle (1986) found both geographic and depth 
stratification of different stages of both phyllosoma and 
pueruli collected during yearly larval tows in the Cali- 
fornia Current Ecosystem, suggesting active ontogenetic 
shifts in larval depth preference. Incorporating such lar- 
val behaviours into biophysical models of dispersal in a 
congeneric species, Panulirus argus, resulted in a 60% or 



greater decrease in the average distance a larva is pre- 
dicted to settle from its release site compared with sim- 
ulations of larvae that remain on the surface (Butler 
et al. 2011). Lobsters in the California Current Ecosys- 
tem rely on kelp forest habitat for survival: P. interrup- 
tus has much higher relative survival rates in the 
presence of kelp than in surrounding areas where kelp 
is absent (Mai & Hovel 2007). Given the ephemeral 
dynamics of kelp forest habitat (Reed et al. 2006) and 
the high variability in ocean conditions in this region, 
these lobsters may have evolved similarly complex 
behaviours as their congeners to increase successful 
local recruitment despite their extremely long PLD 
(e.g. Shanks & Eckert 2005). 

Conclusion 

Here, we present a novel approach to understand 
contemporary drivers of population differentiation in 
systems with high gene flow. In isolation, the population- 
level data present a commonly documented scenario 
among marine species: no evidence for any particular 
regional separation or isolation-by-distance patterns, 
but low and significant pairwise differences among 
populations. While the agreement between nuclear and 
mitochondrial markers confirms that the results are not 
due to statistical artifact, the nature of F-statistics leaves 
us without a clear indication of what is driving the pat- 
tern of genetic differentiation. The addition of kinship 
analyses reveals how alleles are shared between indi- 
viduals, rather than just among populations, and pro- 
vides an independent test of the hypothesis that 
population genetic structure as measured by F-statistics 
is a result of population connectivity. In this case, the 
majority of locations contained an excess of closely 
related individuals. This supports the inference that 
either self-recruitment or some form of coordinated lar- 
val delivery is driving population-level genetic differ- 
ences in a species that would be expected to be broadly 
dispersive throughout its range, given its extremely 
long PLD. In combination with regional oceanographic 
data and larval dispersal behaviour, kinship analyses 
provide evidence for a mechanism of differentiation in 
an otherwise murky population genetic data set. The 
ability to directly test hypotheses about what drives 
population genetic substructure in high gene-flow 
species by independent means, such as kinship or 
coalescent analyses, provides greater confidence in the 
underlying causes of population substructure than sum- 
mary statistics alone. As the ease of developing greater 
numbers of genetic markers increases, individual-based 
analyses such as relative kinship indices can provide a 
valuable complement for understanding patterns in 
traditional population genetics data sets. 
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types and an mtDNA COI sequence alignment: Dryad 
entry doi: 10.5061 /dryad.j6332. 

Supporting information 
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Table SI Isolation by distance (IBD) results from Mantel tests 
conducted in genodive 2.0b20 to examine patterns of genetic 
differentiation across multiple spatial scales based on sampling 
locations: all sites, only California sites, only Mexico sites, only 
continental sites, and only island sites. 

Table S2 Locus-specific diversity indices and F sx values as cal- 
culated in genodive for each of the seven nuclear microsatellite 



loci, and null allele frequency for each microsatellite locus as 
calculated in ML-Relate. 

Table S3 Comparison of ML-Relate calculations of relatedness 
(r) and genodive calculations of kinship for each sampling loca- 
tion. 

Fig. SI Linear regression of mean pairwise local F ST (squares, 
dashed line) and Dest Chao (circles, solid line) at each site on 
the combined proportion of full- and half-sibs as determined 
in ML-Relate for seven microsatellite loci. 

Fig. S2 Log-linear regression of mean relatedness ()') at each 
site on the distance (km) to the nearest edge of an area of high 
upwelling intensity (from Fig. 6a). 

Fig. S3 Median-joining network for Panulirus interruptus 
mtDNA, constructed using 454 base pairs of cytochrome c oxi- 
dase subunit I (COI) from each of 931 individuals in the pro- 
gram NETWORK. 
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